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ABSTRACT 

We report on synoptic observations at 3.6 and 4.5/iin of young stellar objects in IC 348 with 38 
epochs covering 40 days. We find that among the detected cluster members, 338 at [3.6] and 269 
at both [3.6] and [4.5], many are variable on daily to weekly timescales with typical fluctuations of 
^0.1 mag. The fraction of variables ranges from 20% for the diskless pre- main sequence stars to 60% 
for the stars still surrounded by infalling envelopes. We also find that stars in the exposed cluster 
core are less variable than the stars in the dense, slightly younger, south-western ridge. This trend 
persists even after accounting for the underlying correlation with infrared SED type, suggesting that 
the change in variable fraction is not simply a reflection of the change in relative fraction of class I 
vs. class II sources across the cloud, but instead reflects a change in variability with age. We also 
see a strong correlation between infrared variability and X-ray luminosity among the class II sources. 
The observed variability most likely reflects large changes in the structure of the inner wall located 
at the dust sublimation radius. We explore the possibility that these structural perturbations could 
be caused by a hot spot on the star heating dust above the sublimation temperature, causing it 
to evaporate rapidly, and increasing the inner radius for a portion of the disk. Under a number of 
simplifying assumptions we show that this model can reproduce the size and timescale of the 3.6 and 
4.5/im fluctuations. Regardless of its source, the infrared variability indicates that the inner disk is 
not a slowly evolving entity, but instead is a bubbling, warped, dented mass of gas and dust whose 
global size and shape fluctuate in a matter of days. 



1. INTRODUCTION 

Variability has been a defining characterist ic of youn g 
pre-main sequence stars for over 60 years 
These fluctuations occur not only over a wide range of 
timescales, but over a wide range of wavelengths and 
have been used to deduce various properties of these 
systems. Daily to weekly optical fluctuation s up to 1-2 
magnitudes are common (jHerbst et al.l 11994'). Such be- 
havior has been used to infer the presence of hot and cold 
spots covering a small to moderate fraction of the star 
and is often used to stu dy the angular momentum evolu- 
tion o f young stars (e.g. iHerbst "et~aL|[200l iRebull et al.l 
2006) . X-ray variability, often associated with flares from 
magnetic rec onnection events , is also common a mong 
young stars (jWolk et al.ll2005t IStelzer et al.l [2007h . Se- 
lect young stars have even been o bserved to be variable at 
wave l engths as long as the radio (Fcigclso n fc Montmerld 
1999; iForbich et afl[2008l . 

Infrared wavelengths are certainly no exception when 
it com es to young stellar variability. iCarpenter et all 
(2001) find that many stars toward the star-forming 
cloud Orion A are variable in the near-infrared on 
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timescales o f days to weeks. iMorales- Caldcro n et al.l 
(|2009l 1201 If) find that ~70% of youn g stars are variable 
at 3.6 and 4.5/xm and iBarsonv et all (|2005f ) even find a 
significant variable fraction at 10/im on similarly short 
timescales. These wavelengths are particularly interest- 
ing because they probe the terres trial planet forming re - 
gion within a few AU of the star (|D'Alessio et al.i r2006') . 
The variability has the potential to provide information 
about a region of the disk that is well below the resolu- 
tion limit of current telescopes. 

These surveys have demonstrated the ubiquity of in- 
frared variability, while more targeted campaigns of a 
small handful of objects have probed the details of disk 
structure. They have been able to infer vertical pertur- 
bations in the dust at t he very inner edge of the disk 
([Hutchinson et al.l 119941: lJuhasz et al.l 120071 : iSitko et al.1 
2008; Espailla t et al.1 120111 ). These perturbations pro- 
duce a distinctive 'seesaw' behavior in the spectral en- 
ergy distribution (SED) in which the short wavelength 
flux (A < 8//m) increases as the inner disk grows and 
the long wavelength flux decreases due to the additional 
shadowing of the outer disk. It is difficult to reproduce 
a SED that appears to pivot around a point by the ro- 
tation of a vertical perturbation, without also changing 
the n ature of that perturbation (Flah erty fc: Muzerolld 
l2010f) . suggesting that these structures are constructed 
and destroyed on short timescales. Inhomogeneous 
disk struct ure has also been s een in resolved observa- 
tions (e.g. lArnold et al.l 120121 ) and has been inferred 
from rapid occupation events by dust clouds above the 
disk (iBouyier et all 120071: iMorales-Calderon et al.ll2011t 
ISitko et all 120081 1 Alencar et al.ll2010fl . The 'seesaw' be- 
havior appears to be common among disks with large ra- 
dial gaps and holes (Espaill at et al.1 12011: Fla herty et al.l 
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120111 l2012f) but may not be common among less evolved 
objects ( Kospal et al.l [20121 Megeath et al. accepted). 
Among binary systems, variabilit y in the optically 
thin emission ha s been observed ( Skemer et alJ 120101 : 
iNagel et alJ[20ll and can be connected to the streamers 
of dust flowing through the cl eared out region of the disk 
surrounding the binary stars ()Nagel et al.ll2012T ). 

Throughout these targeted studies, infrared variability 
has been used to infer new aspects of the disk struc- 
ture, but only for a limited number of objects. Here we 
present an infrared monitoring campaign of the young 
(~ 3Myr), nearby (~320pc) cluster IC348 in an attempt 
to characterize variability for many young stellar objects 
over a wide range of evolutionary stages. In section 2 we 
describe our Spitzer warm-mission monitoring campaign 
and in section 3 we discuss our technique for selecting, 
and our sensitivity to, variable stars. Section 4 discusses 
variability among the cluster members as well as how 
this variability depends on various system parameters. 
We conclude in section 5 with a discussion of potential 
causes for these fluctuations. 

2. DATA 

We performed repeated IRAC 3.6 and 4.5 jxm observa- 
tions (PID 60160) of IC 348 over the course of 40 days 
from 3 Oct to 7 Nov 2009 to study the infrared variability 
of young stellar objects within this cluster. The cadence 
of these observations was varied from once every four 
hours to once every two days to capture a wide range of 
possible timescales for the variability. We mapped the 
cluster with a 5 by 5 grid, with each grid point sepa- 
rated by 260" and the total field centered at 03:44:20 
+32:03:01. Since the fields of view of the [3.6] and [4.5] 
arrays do not completely overlap, our final map covered 
by both bands is not a square and has a total area of 
roughly 0.4° by 0.25°. We chose a field center that 
is slightly south of the center of the main cluster core 
to cover the YSOs in the more highly obscured south- 
western ridge. The incomplete overlap between the [3.6] 
and [4.5] fields of view, and our choice of field center, re- 
sults in many cluster members only having [3.6] data with 
no corresponding [4.5] data. Images were taken in HDR 
mode with a frame time of 12 seconds and 3 cycles at each 
position (to allow for removal of cosmic rays). The depth 
allowed us to detect sources down to [3.6] , [4.5] ~16 in each 
epoch , sim ilar to earlier Spitze r studies by ILada et alJ 
(2006) and lMuench"eTaLl (120071). The IRAC data re duc- 
tion pipeline is described in] Gutcrm uth et al.l (|2009f ) and 
IMorales-Calderon et al] (j2012ft with updates appropriate 
for the warm Spitzer mission described in Gutermuth et 
al in prep. This monitoring was part of a larger multi- 
wavelength campaign that includes MIPS [24] photome- 
try (PID 40372), IRAC 3-8^m cold mission photometry 
(PID 50596) (Muzerolle et al. in prep), Herschel PACS 
photometry (Balog et al. in prep) as well as ground- 
based spectroscopy (jFlahertv et al.ll20TTl 12012ft . 

3. DEFINING VARIABLES 

When defining a star as variable we strive for a method 
that maximizes sensitivity while still being applicable to 
many stars. We must first accurately characterize the un- 
certainty in the data. We start by calculating the rms, 
(Jobs, for every star that was detected in both [3.6] and 
[4.5] in all 38 epochs with a photometric uncertainty less 



than 0.1. The photometric uncertainty, as reported by 
the source extraction routines, is a combination of shot 
noise, uncertainty in background estimation and mea- 
sured standard deviation in the sky annulus. The actual 
dispersion in the data may be larger than this due to 
detector effects that have not been fully corrected. The 
size of these systematic errors can be measured using the 
rms, (J bs, of the non- variable sources. The a bs is defined 
as: 

2 nY^ =1 Wi{mi — m) 2 
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where n is the number of observations, rrii is the mag- 
nitude on a given epoch, fh is the average magnitude, 
Wi = 1 /erf and a% is the photometric noise. In Figure Q] 
we show the a b s in each band as a function of [3.6] and 
[4.5] magnitude. The main locus contains non-variable 
sources and traces the uncertainty in the data. The un- 
certainty is relatively flat for stars brighter than 13th 
magnitude, but increases sharply for fainter sources. We 
find that the a bs is larger than the photometric uncer- 
tainties (shown as the dashed line) and the difference be- 
tween (Jobs and the photometric uncertainties increases 
from 13th to 10th magnitude. To more accurately esti- 
mate the uncertainties, we fit a line to the locus of non- 
variable sources (shown as a grey line in Figure [T]). For 
each star in our sample, we use its average magnitude to 
interpolate along this line and derive the uncertainty that 
will replace its photometric uncertainty Brighter than 
10th magnitude there are not enough stars to create a 
well-defined locus and we assume the uncertainty is 0.017 
and 0.014 mag for [3.6] and [4.5] respectively. A possible 
concern with this approach is that the nebulous emis- 
sion associated with the molecular cloud creates a highly 
variable background that may lead to non-uniform un- 
certainties across the field. We examine the rms of those 
stars seen against the nebulosity, versus those seen off 
the nebulosity and find no difference in the location of 
the non-variable locus. This indicates that the nebulos- 
ity does not substantially effect our uncertainties and in 
turn does not change our sensitivity to variability. 

With an accurate estimate of the uncertainty in the 
data, we can determine if the spread in the flux of a 
particular star is larger than the uncertainties. This is 
quantified as the reduced chi-squared, xl, which is de- 
fined in each band as: 
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where v = n— 1 is the number of degrees of freedom. Fig- 
ure[2]shows the distribution of xl in each band along with 
the expected chi-squared distribution given the number 
of degrees of freedom. The non- variable sources have xl 
centered at 1 with a spread that matches the expected 
distribution. The long tail toward high xl arises from 
the variable stars and any source with xl > 3 was taken 
to be variable. This boundary is confirmed by visual 
inspection of light curves that fit this criteria. 

The use of the x 2 could be biased by a single data 
point in one band with anomalously high or low flux. 
In fact we find 13 cluster members (LRLL 39, 110, 154, 
186, 286, 298, 329, 353, 478, 1477, 1707, 1719, 1923) for 
which the xl changes by more than a factor of two in 



3 




Figure 1. Observed RMS versus magnitude for all sources, not just cluster members, detected in every epoch for [3.6] (left) and [4.5] 
(right). The photometric uncertainties (dashed line) systematically underestimate the true uncertainty as measured from the non- variable 
stars (solid line) . We use the measured uncertainty of the non- variable stars as an estimate of the true uncertainty in the data. For clarity 
we have excluded the one source with a t, s greater than 0.3 from these plots. 
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residual (assumed to be 1 for our two band data), and 
Pi = Su 1)62(1) is the product of the normalized residuals. 
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The normalized residual in each band is defined as: 
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Figure 2. Reduced x 2 distribution for [3.6] (left) and [4.5] (right). Solid grey line shows the expected distribution if there were no variable 
sources. The substantial population with xt > 3 are all variable in these bands. For clarity we have excluded the 40 and 29 points with 
X 2 , ([3. 6], [4. 5] respectively) greater than 10. 

either [3.6] or [4.5] due to a discrepant data point result- 
ing in a misclassification as a variable star. For these 
sources, as well as the rest of the cluster members, Ta- 
ble [5] lists the membership information (RA, Dec, T e //, 
L», etc.) while Tables I3T41 list the photometry. All show 
either one or two data points, in only one band, that are 
either much brighter or fainter than the rest of the light 
curve. Of these six are class III sources, three are evolved 
disks and four are class II sources (the definition of class 
Ill/evolved/class II is discussed below). The presence of 
this small number of dubious variables will not strongly 
affect our results, and we do not explicitly exclude them 
from our sample. 

Another stat istical tool fo r selecting variables is the 
Stetson index (|Stetsonlll996f) . which takes advantage of 
both bands to look for correlated variability, and has 
been successfull y used in the past to se lect variable young 
stellar objects ([Carpenter et al.ll2001[ ). The Stetson in- 
dex is defined as: 
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where n is the number of epochs for which both [3.6] 
and [4.5] exist, gt is the weight given to each normalized 



Since many of the variable stars show correlated fluctu- 
ations in both bands (see discussion below) this index 
is particularly effective at picking out weak fluctuations 
that would not have been flagged based on the xt- Fig- 
ure [3] illustrates the distribution of Stetson indices, which 
shows a clear locus for the non-variable sources and a 
long tail from the variable stars. Stars with S > 0.45 
are taken to be variable. This boundary is defined based 
on visual inspection of light curves, and confirmed with 
the xt value, which independently indicates that many 
of these stars are variable. 

To search for periodic variables we de trend the ligh t 
curve, followed by analysis with SigSpec (Reegcn 2007). 
We first remove slow trends in the data by subtracting a 
second-order polynomial fit from the light curve. This is 
done to ensure that the mean of the light curve is con- 
stant over the visibility window, which is one of the re- 
quired assumptions of many forms of time-series analysis 
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Figure 3. Stetson index as a function of [3.6] for all sources de- 
tected in every epoch. The dashed line is our adopted boundary 
for splitting non-variable and variable sources. For clarity we have 
excluded the nine sources with a Stetson index greater than 5. 

(jBrockwell fc Davisl I2002T ) . After subtracting the poly- 
nomial fit, w e look fo r significant periods using the rou- 
tine SigSpec (|Reegenll2007l) . Unlike a Lomb-Scargle Peri- 
odogram, SigSpec takes into account both the amplitude 
and phase of the fourier transform of time-series data. It 
analyzes the significance of the amplitude and phase and 
derives a spectral significance, which corresponds to the 
probability that the signal is due to random noise. When 
searching for periods, we consider a significance threshold 
of 5 (a 1 in 10 5 chance that the signal is due to random 
noise) for a signal in one band or a threshold of 4 if the 
signal is present in both bands. Based on simulated light 
curves we find that we cannot accurately recover periods 
longer than 25 days (see discussion below) and we ignore 
any star that has a significant periodic signal longer than 
this limit. 

3.1. Sensitivity 

Given our wide range of fluxes, it is important to un- 
derstand our sensitivity to fluctuations of different sizes 
as a function of magnitude. We can estimate our sensi- 
tivity by simulating light curves with a given fluctuation, 
sampling these light curves at our cadence, adding noise, 
and then calculating the xt and Stetson index to deter- 
mine if they show variations. For the periodic stars we 
start with a sine curve with a period of 15 days with 
[3.6],[4.5]=13, fainter than the majority of cluster mem- 
bers, and the same fluctuation in each band. We create 
10000 stars, add noise to each one, and measure how 
many are detected as periodic as a function of the size of 
the fluctuation. In figure H] we show that we can detect 
all of the stars with large fluctuations, but this sensi- 
tivity severely drops off at small fluctuations. Including 
data from both bands increases our sensitivity by a fac- 
tor of 2, which makes up for the fewer sources that are 
detected in both bands. We estimate that we are able 
to detect >99% of periodic variables with fluctuations 
larger than 0.035 mag. Beyond [3. 6], [4. 5] =14. 5 our sen- 
sitivity severely decreases (Figure [5]) . 

We are not uniformly sensitive to periodicity on all 
timescales because of our limited temporal coverage. 
Rapid fluctuations will not be detected as periodic be- 



cause of our cadence and our finite observing window 
prevents us from deriving long periods. We estimate our 
sensitivity by calculating the period for our simulated 
stars and comparing this derived period to the model pe- 
riod. We find that the stars are not marked as periodic 
if their period is less than 2 days and we underestimate 
the period when it is longer than 25 days (Figure [5]) . For 
this reason we only mark a star as periodic if its period 
lies between 2 and 25 days. 

Estimating the sensitivity to aperiodic variations is 
more complicated. For periodic fluctuations we can use 
a sine curve as the basic model for the light curve, but 
no comparable model has been shown to model the light 
curves of aperiodic fluctuations accurately in young stel- 
lar objects. One potential model is the damped random 
walk model from iKellv et al . (2009) that has been used 
to successfully model quasar variability. This continuous 
first order autoregressive process consists of a white noise 
perturbation that is exponentially damped at a charac- 
teristic timescale. This process is described by the fol- 
lowing stochastic differential equation: 

dX(t) = —X(t)dt + aVdle(t) + bdt (5) 

T 

where r is the relaxation time of the process X(t) and 
e(t) is a white noise process with zero mean and vari- 
ance equal to 1. While this model has been connected 
to accretion disk physics, which may not be applicable 
here given the differences in accretion disks surrounding 
super-massive black holes and pre-main sequence stars, 
it can be used independent of a physical interpretation. 
We employ it simply as a statistical model and do not 
attempt to interpret the various terms as directly relat- 
ing to disk physics. In Figure [7] we show observed light 
curves for three cluster membe rs along with simu lated 
light curves using the model of IKellv et all l|2009[) . We 
are able to generate model light curves that accurately 
match the observations, allowing us to use this model 
to test our sensitivities. While this model has been up- 
dated to include fluctuations on a range of timescales 
(jKellv et alJl201lT >. for simplicity we restrict ourselves to 
the single-timescale model. We construct light curves 
with t = 15 days and then sample those light curves 
with our cadence. Figure 2] shows our sensitivity to vary- 
ing fluctuations for a star with [3.6],[4.5]=13 when using 
just the [3.6] data and when including both data sets. 
As with the periodic fluctuations we find that includ- 
ing both bands increases our sensitivity by a factor of 
2. We estimate that we are able to detect >99% of the 
variables with fluctuations larger than 0.04, similar to 
the periodic variables, although this conclusion is highly 
model dependent. Beyond [3. 6], [4. 5]=14. 5 our sensitivity 
severely decreases (Figure [SJ. 

This model, with its characteristic timescale, allows 
us to test our sensitivity to fluctuations on different 
timescales. Our cadence and total time coverage will 
strongly bias our sensitivity, and we can determine this 
sensitivity by 'observing' fake light curves at the same 
sampling rate as the actual observations. To evaluate 
the fluctuations statistically as a function of timescale 
for the observed light curves, we choose to use the dis- 
crete correlation function, which is similar to the auto- 
correlation function but does n ot assume that the sam - 
pling is evenly spaced in time (|Edelson fc Krolik1 [1988'). 
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Figure 4. Our sensitivity to variables as a function of the size of the variations. We simulate 10,000 13th magnitude stars with a period 
of 15 days (left panel), 5,000 stars with stochastic fluctuations with a timescale of 15 days (middle panel) and 5,000 stars with fluctuations 
with no characteristic timescale (right panel) with a given fluctuation size and measure the fraction of these simulated stars that are 
recovered as being variable using the same criteria as applied to the cluster sample. The solid line shows our sensitivity when both [3.6] 
and [4.5] data are included while the dashed line shows our sensitivity when only using the [3.6] data. Including the second band of data 
increases our sensitivity by about a factor of 2. 

where m,,mj are the magnitudes measured on epochs 
fh is the average magnitude, a is the standard devi- 
ation of the magnitudes, and e is the measurement error. 
By binning over M pairs we create the discrete correla- 
tion function, DCF(r): 
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Figure 5. Sensitivity as a function of magnitude. Simulated stars 
are either given periodic or stochastic fluctuations of 0.04 mag, 
and the fraction recovered is measured as a function of magni- 
tude. Different lines refer to different models of the underlying light 
curve (see text for details). Our sensitivity severely drops off at 
about 14th mag. 
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Figure 6. Comparison of the period derived from the peri- 
odogram versus the given period for fake light curves. Error bars 
show the dispersion in the derived period. Above 25 days the mea- 
sured period underestimates the actual period of the light curve 
and we do not mark stars as periodic if we derive a period longer 
than 25 days. 

To calculate the discrete correlation function we first cal- 
culate the unbinned discrete correlations, defined as: 



UDCFi, 
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DCF(r) = ^UDCF tJ 
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We create simulated light curves as before, but instead 
of simply measuring the size of the fluctuations in these 
light curves we now calculate the DCF for each curve. 
To examine the general shape of the DCF for a particu- 
lar timescale we average the DCF derived from each of 
our 5000 simulated light curves. We show the results 
for different timescales in figure |U The DCF for light 
curves with timescales below ^2 days are completely flat. 
Above 25 days the shape of the DCF does not change, 
indicating that we are not able to constrain timescales 
longer than 25 days. Also in figure [8] we show the DCF 
for individual simulated light curves with a timescale of 
20 days. The noisiness in each individual DCF indicates 
that we will be unable to use the observed DCF for the 
IC 348 cluster members to constrain the timescale accu- 
rately. In figure |9] we show the DCF for a sub-sample 
of cluster members. These data also show the noisiness 
that prevents us from accurately deriving a timescale for 
the variability. In general the DCF are not completely 
flat, indicating that the light curves vary on timescales 
longer than 2 days. Increasing the number of epochs, as 
well as the time baseline, would improve our ability to 
constrain the timescale of aperiodic fluctuations. 

We can begin to probe the fluctuations on longer 
timescales by comparing our fall 2009 warm-mission 
[3.6] and [4.5] photometry with our March 2009 IRAC 
cold- mission da ta (Muzero ll e et al. in prep) as well 
as the GTO (lLada et all I20M PID 6) and C2D 
(|Jorgensen et al II2006L PIP 178) data taken in 2004. We 
find that 96%(92%) of the cluster members have a mean 
[3.6] ([4.5]) magnitude from March 2009 between the min 
and max magnitudes of the fall 2009 data. Only eight 
stars (LRLL 31, 35, 99B, 245, 382, 1683, 9099), 3% of 
the cluster, have both March 2009 [3.6] and [4.5] beyond 
the min and max of what is seen in our fall 2009 data. 
Roughly 6% of the cluster members have both [3.6] and 
[4.5] in the GTO or C2D data outside of the range seen 
in our fall 2009 data. This suggests that large year-long 
fluctuations are uncommon among the cluster members, 
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Figure 7. [4.5] light curves of LRLL 15, 26 a nd 31 (red c r osses , left to right) along with stochastic models (black squares) 
timescale of 20 days. The stochastic model of IKellv et al.l ( [20091 ) is able to accurately reproduce the observed light curves, 
it may be an accurate description of the underlying process that drives the variability, and can be useful for understandinj 
of our survey 
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Figure 8. On the left is the median discrete correlation function (DCF) for 5000 stars with aperiodic variability at different timescales. 
Light curves with timescales between 2 and 25 days have appreciable differences in the DCF. On the right we show the DCF for five light 
curves with a 20 day timescale (colored solid lines) along with the median DCF of 5000 stars with this timescale. The large variability in 
the shape of the DCF between different stars demonstrates the difficulty in deriving the variability timescale for a single star, despite the 
change in average shape of the DCF with timescale seen in the left-hand panel. 



o 

Q 




20 

Log (doys) 



40 



Figure 9. Discrete correlation function for sample young stellar 
objects. While the shape of the DCF depends on the variability 
timescale, the limited sampling of our data makes the observed 
DCF too noisy to draw strong conclusions about the timescale. 

although more detailed analysis is needed. 

Using the damped random walk model to estimate our 
sensitivity to fluctuations of different timescales and sizes 
is effective, but gives us a result that is highly model 
dependent. We can instead simulate the stochastic vari- 
ables using a different model that has no characteristic 
time scale. The stochastic model described above has a 
power-spectrum density (PSD) that is a power law over 
many fr equencies, but has a break at a characteristic fre- 
quency (|Kellv et al.ll2009f) . We can also model the PSD 



with a power law (w~ 2 ) with no break (jTimmer et al.l 
f 995). This model also produces light curves that reason- 
ably match our data and can be used to estimate the sen- 
sitivity. Light curve simulations are performed as before 
and the results are shown in Figure |4j We estimate that 
we are able to detect >99% of the variables with fluctu- 
ations larger than 0.045, similar to the previous model. 
Our sensitivity also drops off beyond [3. 6], [4.5] =14.5, as 
was seen with the stochastic light curve model. 

Overall, we find that we are sensitive to almost all fluc- 
tuations larger than 0.04 mag, with decreasing sensitivity 
beyond [3. 6], [4. 5] =14. We can detect periodic fluctua- 
tions if the period is between 2 and 25 days, with dimin- 
ished sensitivity to periods below 2 days and a systematic 
underestimation of the period when it is above 25 days. 
We have very little sensitivity to the exact timescale of 
aperiodic fluctuations, although there appear to be few 
hourly fluctuations or large year to year fluctuations. 

4. CLUSTER VARIABILITY 

With an understanding of the sensitivities and lim- 
itations of our data , we examine the IC 348 cluster 
(|Luhman et al.lr2003[ ) . The cluster m e mbers h ip listed has 
been co mpiled from Luhm an et al.l (I2003D : ILada et ahl 
(120061); iMuench et alj (|2007| ): IMuzerolle et alj (j2006fk 
ILuhman et alj (|2005D with the membership information 
(LRLL #, RA, Dec, T e// , etc.) listed in Tableland the 
[3.6] and [4.5] photometry listed in Tables [3] and 0] re- 
spectively. We detect 338 of 377 cluster members at [3.6] 
and 269 in both bands. The distribution of fluxes for 
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the detected cluster members is shown in Figure [TU1 The 
majority of the cluster members that are not detected at 
[3.6] are either outside the field of view (24 sources), are 
in a very early evolutionary state or have very low intrin- 
sic luminosity that make them faint at these wavelengths 
(6 members), are close to another star, which confuses 
their photometry (3 members) or have very larger pho- 
tometric uncertainties (> O.lmag, 2 sources). Due to 
the placement of the field, and the incomplete overlap 
between the [3.6] and [4.5] fields of view, there are a sig- 
nificant number of cluster members detected at [3.6] but 
not [4.5]. Only one source, located at the very edge of the 
field of view, is detected at [4.5], but not [3.6], most likely 
due to imperfect overlap between the two fields of view. 
Most of the cluster members are brighter than 14th mag 
in both bands, which is roughly the magnitude at which 
the uncertainties start to increase dramatically. Of the 
cluster members detected at [3.6], 84 are variable and 254 
are not variable, based only on the [3.6] data. For the 
cluster members detected in both bands, 125 are variable 
and 144 are not variable. In Table [2] we mark the cluster 
members that we consider variable and non-variable, as 
well as whether or not this conclusion was derived from 
one band versus using both bands. The large difference 
in variability fraction reflects the large increase in sensi- 
tivity when both bands are examined together. For the 
rest of our analysis we consider only the 269 stars that 
have been detected in both bands because the increase in 
sensitivity makes up for the decrease in the sample size. 
When discussing variability trends, we will often distin- 
guish between periodic variables and irregular variables. 
Keep in mind that when we say 'irregular variables' we 
really mean to say 'stars that do not have fluctuations 
larger than 0.035 mag on a period between 2 and 25 
days.' Some of the 'irregular variables' may in fact be 
periodic, but at a level or on a timescale that we are not 
sensitive to. 

4.1. What influences infrared variability? 

One distinct advantage of studying the IC 348 cluster 
is the wealth of ancillary data that is available in the lit- 
erature. Its compactness and prox imity make it an idea l 
target for optical spe ctroscopy (Luhman et ail [2 003). 
infrar ed photometry (lLada et al.l 120061: iMuench et al.l 
120071) . X-ray imaging (jPreibisch fc Zinneckerll2002ll2004fl 
and optical photometric monitoring (I Cohen et al.1 12004J : 
iCieza fc Balibeifl2006tlNordhagen et al.l2006t ). These ob- 
servations have led to a well-defined membership sam- 
ple that is complete down to brown dwarf masses, along 
with measures of the Ha emission line strength, infrared 
excess and X-ray luminosity for hundreds of members. 
Determining which of these system parameters influence 
the presence and strength of variability can constrain the 
physical model of the variability. 

The parameters we examine, listed in Table [3J are a 
combination of values from the literature and those de- 
rived from our data. Stellar parameters, such as effective 
temperature (T e ff), stellar luminosity (L„), Ha equiva- 
lent width (EW), extinction (Ay), and X-ray luminos- 
ity (Ly) are all taken from th e literature (jLuhman et al.l 
120031: LPreibisch fc Zinneckerll2002D . To characterize the 
shape of the excess we use the slope, a, of the infrared 
spectral energy distribution (SED), defined as XF\ cx A" 
where a is measured over the 3-8/im IRAC bands. The 



IRAC photometry was measured from a stack of the 
five consecutive days of monitoring taken in March 2009 
(Muzerolle et al. in prep). We favor the added sensitiv- 
ity from this deep stack in exchange for the possibility 
that the shape of the SED changes over time, which is 
occasionally but not frequently seen (Muzerolle et al. in 
prep) . Infrared fluxe s have been deredden ed according to 
the extinction law of Flahert y et al.l (|2007l ) for those stars 
with measured extinction. The value of amAC ranges 
from a = —3 for purely photospheric emission (Class 
III objects), —2.56 < a < —1.8 for evo lved disks, and 
a > -1.8 for full disks ([Lada et al.ll2006| ). Full disks can 
be further divided into class II sources (—1.8 < a < 0) 
and class I sources (a > 0). Unless otherwise stated 
the variable ctiRAC refers to the slope of the dereddened, 
rather than observed, SED. 

We can also characterize the excess based on the dif- 
ference between the observed and the photospheric flux 
at 3.6/im. A photosphere template is normalized to the 
dereddened J band flux and subtracted from the dered- 
dened [3.6] flux (F(3.6) o6s /F(J) ?6s -F(3.6) p/lot /F(J) pho4 ) 
to produce a stellar flux normalized 3.6/xm excess. For 
the photosphere template we use a median of the flux 
from stars within IC 348 with no disk (olirac < —2.56) 
and an effective temperature within 5% of the target star. 
We prefer using th e cluster members to define t he photo- 
sphere because the lKenvon fc Hartmannl (|1995l ) observed 
colors, as well as Kurucz model atmospheres, tend to 
overestimate the photospheric flux beyond K-band and 
thus underestimate the 3.6/xm excess. For the hottest 
stars, where there are few stars with similar e f fective 
temperatures, we use the iKenvon fc Hartmannl (1995) 
observed colors. Both Chirac and the [3.6] excess trace 
the strength of the infrared emission, but do so in slightly 
different ways. The [3.6] excess, despite the larger uncer- 
tainties (~ 0.01), more directly traces the same part of 
the disk as the [3.6] and [4.5] photometry. A system with 
a large excess at 8/xm, but no excess at 3.6/rni would have 
a large olirag but the lack of dust capable of emitting 
in the [3.6] band may lead to less variability. 

The last parameter that we add is the position of the 
young stellar object within the cluster. IC 348 is com- 
posed of an exposed cluster core in the north-east and a 
highly obscured region in the south-west where most of 
the younger class I sources, along with al l of the molec- 
ular outflows and sub-mm cores, reside ([Herbstl l200l . 
The presence of these less-evolved sources indicates that 
the south-western ridge is a region of very recent star- 
formation. We characterize the position as the distance 
from HH 211, a prominent Herbig-Haro object in this 



region. 



4.1.1. Irregular Variables 



During most of our analysis we only consider the irreg- 
ular variables. The periodic variables represent a small 
subsamplc of the entire survey (121 irregular variables 
versus only 13 periodic variables, Table [2]), and possibly 
have a different physical origin to their variability, which 
may bias our results. We include them in the figures, 
but consider the plausible causes of their variability sep- 
arately in section 4.1.3 

We start by considering the shape of the dereddened 
infrared SED, uirac and how it correlates with the frac- 
tion of variables (f va r) and the size of the fluctuations in 




Figure 10. Distribution of the [3.6] and [4.5] magnitudes for the cluster members detected on more than one epoch. The solid, unfilled 
histogram shows the stars detected in at least that band, while the dashed, filled histogram shows the stars detected in both bands. There 
are a substantial number of cluster members detected at only [3.6], but only one star was detected at [4.5] and not [3.6]. 



variable stars for the 239 cluster members with measur- 
able variability and ajRAc- We use linear regression to 
determine if there is a statistically significant correla- 
tion between Chirac and infrared variability. If the slope 
of the linear fit is significantly different from zero, then 
the two variables are correlated. This does require that 
we bin the stars in order to calculate i var and convert 
f var into a logit (=ln(f, mr /(l-f ? , nr ))) in order to make 
it am enable to linear regression (jHosmer fc Lemeshowl 
2000). The 239 stars with greater than two epochs of 
both [3.6] and [4.5] photometry, needed for measuring 
variability, and measurable infrared SED slopes, are split 
into four bins, each with an equal number of sources. The 
parameter f var is defined as the fraction of stars that ex- 
ceed our variability criteria (x 2 > 3 in either [3.6] or [4.5] 
or S>0.45) and the size of the bin in parameter space 
is varied in order to ensure that each bin has an equal 
number of cluster members. To quantify the size of the 
variability we sum the xt m each band, which allows us 
to utilize information in both infrared bands. Again we 
perform linear regression to evaluate the strength of the 
trend statistically. 

We find a strong correlation between f var and o.irac 
(Figure lllb ) with the variable fraction increasing from 
20% for those stars without disks to 60% for stars with 
full disks/envelopes. We also show the cumulative frac- 
tion, although it isn't used to look for trends statistically, 
to visualize the correlation without the need for binning 
(Fig lllb). We also find a significant trend between the 
size of the variability and the SED slope (Figure ITTb). 
The xti while useful statistically, is hard to interpret in 
terms of actual changes in the [3.6] and [4.5] magnitude. 
To more intuitively understand the size of the fluctua- 
tions in Figure 1111 1 we calculate the max-min magnitude 
for each irregular variable and show the median of these 
values in each bin. These bins are the same as used when 
defining the variable fraction, with an equal number of 
cluster members, but not an equal number of irregular 
variables, in each bin. The typical size of the fluctua- 
tions increases from ^0.15 mag to 0.3 mag from diskless 
to class I sources. 

This trend is co nsistent with previous studie s of in - 
frared variability (|Morales-Calderon et all 120091 1201 If ). 
Splitting the cluster members into infrared classes us- 
ing the dereddened 3-8/im SED slope, we find that 
79% of class Is (aiRAC > —0.5), 61% of class lis 



(-1.8 < amAc < -0.5), 44% of evolved disks (-2.56 < 
aiRAC < -1.8) and 20% of class Ills {ptjRAC < -2.56) 
are variable (24,101,64,75 total stars in each group). 
The YSOVAR survey of Orion found variable fractions 
of 84% 70% and 45% for class I, II and III respec- 
tively ()Morales-Calder6n et al.l 1201 If) . These variable 
fractions show the same trend with evolutionary stage 
as in IC 348, although the fractions are higher in each 
SED class. Orion is younger than IC 348, and it is pos- 
sible that as stars evolve they become less variable, a 
possibility discussed in more detail below. In IC1396A 
the variable fractions for class I and class II sources 
are 75% and 58% respective ly, similar to our results 
(jMorales-Calderon et aLll2009f ). 

This analysis can be extended to the other system pa- 
rameters, but we must be careful to avoid degeneracies 
with the shape of the infrared SED. For example, the Ha 
line s trength broadly t races the accretion rate onto the 
star (|Natta et al. 2004) and will be correlated with the 
infrared excess. This may result in a statistically signif- 
icant trend between variable fraction and Ha EW that 
may simply be a reflection of the underlying trend with 
infrared SED shape. To avoid this degeneracy we can 
perform linear regression between the logit or xt and a 
linear combination of system parameters. This amounts 
to looking for a trend between Ha EW and i var or xl 
at a given value of aiRAC- Using a linear combination 
as a simple statistical model to relate i var or xt and 
system parameters almost certainly underestimates the 
complex relationship between different parameters but it 
is useful as a first order approximation to determine the 
most important factors in the system. The F-statistic 
is then used to determine the significance of adding an- 
other variable to our statistical model. While we do con- 
sider sources across all SED classes when performing this 
statistical analysis, we only show the class II sources (- 
1.8< aiRAC <0) in the plots below (Figure [1211201) for 
clarity. 

One possibility is that the correlation with Ha EW, 
or some other parameter, is more fundamental than the 
correlation with SED shape and that we should instead 
look for a trend between Chirac and variability after con- 
trolling for the correlation with Ha EW instead of the 
other way around. We have chosen to use the correlation 
with aiRAC as the more fundamental trend because it is 
the strongest correlation in our data, it has been seen by 
other groups, and requiring an infrared detection puts 





Figure 11. Variability as a function of the shape of the infrared SED (239 members), (a) Fraction of cluster members that have irregular 
fluctuations (black squares) or periodic variability (red triangles) . Each bin contains an equal number of cluster members and the horizontal 
error bars represent the size of the bin while the vertical error bars represent the poisson noise on l va r- (b) Cumulative fraction of irregular 
variables (c) Size of the fluctuations, measured using the sum of the xt m each band, for the irregular variables (black crosses, 120 members). 
Squares show the mean fluctuation in bins with equal numbers of cluster members, but not necessarily an equal number of variable stars, 
(d) Mean of the peak to peak fluctuations in magnitude of the irregular variables at [3.6] (black squares) and [4.5] (red stars). Horizontal 
error bars designate the bin size and vertical error bars are the error on the mean. The bin with the smallest ainAC does not have enough 
variable stars to define a mean, and is assigned a value of zero. 



fewer limits on the number of sources we can study than 
requiring e.g. an Ha EW measurement. With the rest 
of the system parameters, we first examine if there is a 
significant correlation with variability, and then examine 
whether or not this trend persists after controlling for 
the underlying correlation with SED shape. 

Interestingly we find a significant trend between the 
position within the cluster and the variable fraction 
(Fig [T2"a.b) for the 249 cluster members with accurate 
positions and measurable variability. As before, i var is 
defined as the fraction of stars with xt > 3 in either 
[3.6] [4.5] or S> 0.45, separating out stars with evidence 
for periodicity, in bins whose boundaries are set so that 
each bin contains the same number of cluster members. 
The variable fraction decreases from 70% near HH 211 
to 20% on the far side of the cluster (Fig [l2"h). This 
trend is still significant (F=9.5) even when including the 
underlying correlation with aiRAC (Figure H^H.e). The 
trend between t var and position is not due to a change 
in the relative fraction of class I/class II/class III be- 
tween different regions of the cluster. There is a change 
in systems of a given SED shape that causes the vari- 
ability to diminish from one side of the cluster to the 
other. Even when examining only the class II sources, 
l var drops from 70% to 50% across the cluster. This can 
be seen in Figure [JJ2 where we show a map of the stars 
with disks {olirac > —1-8) within the cluster. In an- 
nuli that increase from 5 to 15 arcminutes from HH 211 
the variable fraction decreases from 76% down to 52%. 
Over this SED range i var is fairly constant as a function 
of SED shape (Fig ITTk). but there is still a correlation 
with position. There must be some change in the stars 
with disks that causes them to be less variable the farther 
they are from the dense, embedded region in the south- 



east. There is a hint of a trend between the size of the 
variability and the location within the cluster (Fig|T2"b). 
although this is not as significant. 

We also find a significant correlation between the 
strength of the variability and the effective temperature 
of the star (among 235 cluster members) with higher 
T e ff stars showing larger fluctuations (Figure IT4"b.d). 
even when taking into account the underlying trend with 
Chirac (F=7.7, Fig[T4b). Stars earlier than K0 have fluc- 
tuations of ~ 0.3 mag, while stars later than M2 have 
fluctuations of ~ 0.15 mag (Fig [T4h). There is no corre- 
sponding trend between the logit and T e ff (Fig iMk.b). 
It is likely that the observed relation is an observational 
bias related to the fact that the fraction of stellar flux in 
the [3.6] and [4.5] bands increases as the effective tem- 
perature decreases. The increased fraction of stellar flux 
will dilute the signal from the disk, making it appear as 
though the fluctuations are smaller. We can demonstrate 
this effect by considering a stellar SED modeled with 
a single temperature blackbody underneath disk emis- 
sion at Tdust=1500K. The disk is chosen to have a single 
temperature because the disk emission at [3.6] and [4.5] 
is dominated by fl ux from the inner wa ll at the dust 
destruction radius (jMuzerolle et al.l 12003). The fiducial 
disk flux level is set so that r/f, the ratio of disk flux 
to stellar flux at K-band, is 0.4 and we vary the excess 
about this fiducial value for effective temperatures rang- 
ing from 8000 K (A7) down to 2700 K (~M6). In this 
simple model we implicitly assume that inner rim scale 
height does not depend on stellar mass. Figure [TBI shows 
that stars with smaller effective temperature have smaller 
variations in their [3.6] and [4.5] magnitude for a given 
change in the disk flux, consistent with the effect that 
we observe in the data. Because this trend is an observa- 
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Figure 12. (a) Variable fraction as a function of position (249 members), (b) Cumulative fraction, (c) Mean peak to peak fluctuations 

(d) Fraction of irregular variables split into class II sources (black squares, 95 members) and class III sources (blue diamonds, 84 members). 

(e) Cumulative fraction of irregular variables among the class II sources. There is a significant trend between i va r and position even when 
looking at only the class II sources. 
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Figure 13. Map of cluster members with otmAC > —1-8 (ie. 
stars with disks, 95 members) marked as crosses, with variable 
stars marked with diamonds. The position of HH 211 is marked 
with a red plus sign surrounded by circles at a distance of 5, 10 and 
15 arcminutes. The fraction of variable stars with disks is noted 
at each radius. There is a significant decrease in disk fraction 
with radius even at a fixed SEE) shape, indicating that variability 
depends strongly on position within the cluster. 

tional bias, we do not account for it while searching for 
a physical cause of the variability. 

For the other parameters that we consider (L*, Ha 
EW, [3.6] excess, Ay) we find no statistically significant 
evidence for a correlation with the logit or the x 2 , es- 
pecially after taking into account the underlying trend 
with ai R Ac (Figures [T6)l7ll8|19p . We find that, besides 
the trend between T e ff and size of the fluctuations, vari- 
ability does not depend strongly on the base stellar prop- 
erties of T e f f and L» . This suggests that there is little 
dependence of infrared variability on stellar mass from 
GO to M6 (~ 2.5 - 0.1M o ) The large scatter in f var and 
xt still leave some room for small trends with these pa- 



rameters, but we can rule out large trends (Af var > 20%, 
A[3.6], [4.5] > 0.2). The main limitation in our ability to 
measure the dependence of variability on these parame- 
ters is due to uncertainties in their measurements, which 
can smear out any small trends. It may also be possible 
that observational biases, such as the inability to detect 
Ha in deeply embedded class I sources, restrict our abil- 
ity to detect any real trends in the data. For the rest of 
our analysis we focus on the trends strong enough for us 
to detect, but note that we cannot completely rule out 
correlations with the other parameters. 

One potential bias in our analysis is the fact that the 
source of the variability is likely different between class II 
and class III sources. The variability in class II sources is 
most likely related to the structure of the disk and there- 
fore might be correlated with different system parameters 
than the variability of the class III sources, which is most 
likely driven by star spots. Our method of regression be- 
tween the logit or xt and a linear combination of system 
parameters includes cluster members of all SED shapes 
and a lack of a trend among class III sources may dilute a 
real trend among class II sources. To eliminate this bias 
we repeat the analysis performed above but only looking 
at the class II sources. We use the same definitions of 
i var and variability strength as before, we simply change 
the subset of stars that are examined. We find many of 
the same correlations as before; Chirac is correlated with 
both i var and the variability size, position is still strongly 
correlated with i var and T e / / is still strongly correlated 
with the size of the fluctuations. Interestingly we find 
that f var is marginally correlated with X-ray luminosity 
and the strength of the variability is strongly correlated 
with X-ray luminosity (Figure l20e.f.g). This was not 
the case when including the class III sources (Fig[20k.b). 
indicating that X-ray luminosity is only related to the 
variability of the disk. As the X-ray luminosity increases 
from 10 28 to 10 30 erg cm~ 2 s _1 the variable fraction in- 
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Figure 14. (a) Variable fraction as a function of T e // (235 members), (b) Cumulative fraction, (c) Mean peak to peak fluctuations for 
irregular variables . (d) Size of the fluctuations for all irregular variables (86 members), (e) Size of the fluctuations for only the class II 
sources (86 members). There is a significant trend betweenxj and T e ff even when looking at only the class II sources, but this is most 
likely an observational bias. 
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Figure 15. The variability in [3.6] (left) and [4.5] (right) for varying levels of discrepancy away from a fiducial excess for stars at various 
effective temperatures. As the effective temperature of the star decreases the size of the fluctuations decreases for a given change in the 
excess. This demonstrates that the observed correlation between variability strength and effective temperature is not a physical effect, and 
is simply due to the different relative strengths of the stellar and disk flux as the peak of the stellar SED moves to redder wavelengths. 



creases from 40% to 100% and the size of the fluctuations 
increases from 0.15 to 0.4 mag. For the class III sources 
i var is roughly constant with X-ray luminosity. The same 
trend can be seen us i ng X-r ay fluxes from the recent anal- 
ysis of iStelzer et al.l (|201 2f ) who present a stack of X-ray 
data from multiple Chandra observations that is able to 
detect an additional 30 cluster members. 

4.1.2. Other Trends 

The extensive optical monitoring of this cluster allow 
us to compare the presence of optical variability, in par- 
ticular irregular fluctuations, with the infrared variabil- 
ity. Of the 16 stars marke d as irregular optical vari- 
ables by iCohen et al.l (|2004f ) , all of which have an in- 
frared excess, 88% are variable in the infrared, which 
is higher than when considering the entire sample. Of 
the stars that do not show significant irregular optical 



variability, 33% of the strong disks are variable, 30% of 
the evolved disks are variable and 21% of the non-excess 
sources are variable in the infrared. For the strong disks 
the variable fraction is significantly smaller than the en- 
tire sample(~ 60%) indicating that the irregular optical 
variability and the infrared variability are strongly re- 
lated. 

We also find that the fluctuations at [3.6] and [4.5] 
are highly correlated. This is not surprising given that 
the disk flux at both wavelengths is dominated by emis- 
sion from the puffed up wall at the dust sublimation ra- 
dius. To quantify any variations in color we calculate the 
slope of the data in a plot of [3.6] vs. [4.5], taking into 
account the uncertainty in the measurements and any 
intrinsic sca tter in the relation between [3.6] and [4.5] 
(|Kellvlr2007|) . Slopes greater than 1 corresponds to the 
star getting bluer as it gets fainter (or redder as it gets 
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Figure 16. (a) Variable fraction as a function of L, (223 members), (b) Cumulative fraction of irregular variables, (c) Mean peak to 
peak fluctuations among irregular variables (99 members), (d) Size of the fluctuations. No significant trend is seen between L„ and the 
measures of variability (the size trend goes away when including T e ff). 
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Figure 17. (a) Variable fraction as a function of Ha EW in units of A (154 members), (b) Cumulative fraction, (c) Mean peak to peak 
fluctuations for the irregular variables, (d) Size of the fluctuations (72 members), (e) Variable fraction split into class II (59 members) 
and class III (56 members)sources. (f) Cumulative fraction of irregular class II variables, (g) Size of the fluctuations among class II stars. 
While there is an increase in i var with Ha emission, this trend disappears when correcting for the trend with Chirac- 
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Figure 18. (a) Variable fraction as a function of [3.6] excess (F(3.6) i, s /F(J) (, s -F(3.6) p h i/F(J)phoi) (126 members), (b) Cumulative 
fraction, (c) Mean peak to peak fluctuations for the irregular variables, (d) Size of the fluctuations (70 members), (e) Variable fraction 
split into class II (63 members) and class III (30 members) sources, (f) Cumulative fraction of irregular class II variables, (g) Size of 
the fluctuations among class II stars. While there is an increase in f va r and xi with [3.6] excess, this trend is severely diminished when 
correcting for the trend with amAC- 




Figure 19. (a) Variable fraction as a function of Ay (222 members), (b) Cumulative fraction, (c) Mean peak to peak fluctuations for 
the irregular variables, (d) Size of the fluctuations (98 members). There is no strong trend between Ay of a star and infrared variability. 
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Figure 20. (a) Variable fraction as a function of X-ray luminosity in units of 10 28 erg s _1 (154 members), (b) Cumulative fraction, (c) 
Mean peak to peak fluctuations for the irregular variables, (d) Size of the fluctuations (58 members), (e) Variable fraction split into class 
II (33 members) and class III sources (60 members), (f) Cumulative fraction of irregular class II variables, (g) Size of the fluctuations 
among class II stars. When looking at all of the cluster members together there is no trend between { va r and X-ray luminosity (first two 
panels) but when restricting the analysis to just the class II sources (middle panels) there is a strong trend between variability and X-ray 
luminosity. 

brighter) while a slope of less that 1 corresponds to the 
system getting redder as it gets fainter (or bluer as it 
get brighter). Figure [21] shows the slopes for the star 
with the most highly correlated fluxes. Of these stars, 
76% have a slope within 3<r of 1, indicating no change in 
color, while the rest exhibit significant changes in color. 
Four sources (LRLL 21, 31, 67, 90; see Table [2] for posi- 
tions and other ancillary information) get redder as they 
get brighter, with slopes larger than 1. This is most likely 
due to the larger stellar flux in [3.6] versus [4.5] diluting 
some of the change in [3.6] due to the change in disk 
flux when compared to the change in [4.5]. The eight 
sources that get bluer as they get brighter (LRLL 13, 
75, 245, 276, 435, 1889, 40182, 54361) all have a strong 
infrared excess {olirac > —1-0) and possibly have an en- 
velope contributing scattered light to the [3.6] and [4.5] 
bands. If the variability consists of a short pulse, then 
as the light echo passes through the envelope the level 
of scattered light will increase, as is seen in LRLL 54361 
(Muzerolle et al. submitted). 
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Figure 21. Slope of [4.5] vs [3.6] for those stars where the uncer- 
tainty in the line fit to [4.5] vs [3.6] is less than 20%. Most stars 
show no change in color, although a small handful either get redder 
as they get fainter or get bluer as they get fainter. 



4.2. Periodic Stars 

We detect 13 stars that show periodic behavior in their 
infrared variability, based on the SigSpec analysis of ei- 
ther [3.6] or [4.5] (Table [T]). Another 14 stars have weak 



evidence for periodicity based on SigSpec but are likely 
periodic because either the same period is derived in both 
[3.6] and [4.5] or the infrared period matches the optical 
period. That is, about 25% of the variable stars are 
likely to be periodic; the fluctuations in the majority of 
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cases are dominated by short-term aperiodic variations. 
When the [3.6] and [4.5] period do not agree we con- 
sider the period with the strongest significance as the 
actual period. Of the 13 certain periodic stars, the pe- 
riods range from 2 to 22 days, as expected given our 
sensitivities. Eight of the stars are periodic in the op- 
tical, and of these five (LRLL 64, 69, 87, 116, 125; see 
Table [2] for positions and other ancillary information) 
have infrared periods that match their optical periods 
suggesting that the same physical process, most likely a 
cool spot rotating across the star, is driving the infrared 
and optical variability. All of these stars have a weak 
excess {aiRAC <-2.3), consistent with the [3.6] and [4.5] 
fluxes being dominated by photospheric emission as well 
as fluctuations small enough (er < .050) to be consis- 
tent with cool spots. Of the three stars whose optical 
and infrared periods differ, LRLL 29 has an infrared pe- 
riod less than the optical period while LRLL 36 and 90 
and infrared periods much longer than the optical period. 
LRLL 29 is a class III source (chrac — —2.73) with a 
small fluctuation (a = 0.03) and we are likely seeing a 
cool spot at a different rotation period than previously 
measured. For these last two sources, both have strong 
infrared excesses that would mask out stellar fluctuations 
and the variability is most likely due to a long-lived struc- 
ture within the disk. The periods (14.7 days for LRLL 
36 and 20 days for LRLL 90) corresponds to material at 
^0.1 AU from the central star. Of the other sources that 
do not show strong periodic fluctuations in the optical, 
LRLL 94 does not have an excess, and we are likely see- 
ing rotating spots, while the other three sources (LRLL 
15, 186, 435) all have a strong excess and the periodicity 
is likely due to rotating structure within the disk. The 
star LRLL 54361 is notable not only because it is a pe- 
riodic class I source, but also because it is the strongest 
variable in our sample, with fluctuations up to 2.5 mag- 
nitudes in both [3.6] and [4.5]. We refer to Muzerolle et 
al. submitted for more details on the observations, and 
possible physical cause, of the variability in this source. 

4.3. 'Dippers' 

IMorales-Calderon et al.l (|2011h identified a subsample 
of infrared variables in Orion that exhibit brief dips in 
their light curves. These flux drops of a few tenths of a 
magnitude are seen in both [3.6] and [4.5] as well as con- 
temporaneous optical and near-infrared data, and last 
for typically 1-2 days. In some cases they are periodic, 
recurring on weekly timescales, but sometimes only one 
event is seen. The fact that the dips are more extreme at 
shorter wavelengths indicates that they are most likely 
due to extinction events in which the star is temporarily 
obscured by a clump of dust. The rapid timescale in- 
dicates that the clump of dust resides in the disk, close 
to the inner edge. These exti nction events are si milar 
to those seen around AA Tau (IBouvier et al.l 12007ft and 
the c lass of object known as UXors ((Waters fc Waelkensl 
[1991 . 

We do not see any sources that we can confidently clas- 
sify as dippers in our sample, but this may be due to 
our sparse observing cadence. We can simulate our sen- 
sitivity to dippers by randomly generating light curves 
with dips whose chara c terist ics match those seen by 
IMorales-Calderon et al.1 (|201lD and running these arti- 
ficial light curves through our algorithms for detecting 



variables. Our simulated cluster has the same [3.6] and 
[4.5] magnitude distribution as the actual observations in 
order to match our measured sensitivities. We randomly 
select stars to be 'dippers' and assign them a dip with a 
strength between 0.05 and 0.2 magnitudes and a dura- 
tion between 0.5 and 2 days. We repeat this simulation 
1000 times for a given 'dipper' fraction and determine 
how many of these artificial clusters have more than one 
detected dipper. Our 99% confidence limit corresponds 
to the dipper fraction for which fewer than 1% of the 
artificial clusters have more than 1 dipper. In this way 
we can put an upper limit of 5% on the dipper fraction. 
This is roughly consistent wit h the fraction derived by 
IMorales-Calderon et all (|201lD . 

5. WHAT IS HAPPENING TO THE DISK TO CAUSE THIS 
VARIABILITY? 

Among the IC 348 cluster members we find that many 
of the stars show infrared variability. This variability 
is more common among stars with a strong infrared ex- 
cess, as well as stars with a large X-ray luminosity and 
stars located in the south-western ridge. Stars with ir- 
regular optical fluctuations are much more likely to be 
infrared variables than those without such fluctuations. 
Very few of the variables are periodic, or appear to have 
flux changes consistent with rapid extinction events. The 
[3.6] and [4.5] flux is highly correlated for the majority 
of stars, suggesting a common origin for their variability. 
We now discuss what could be causing this variability. 

5.1. Variable Structure 

Our observations at 3.6 and 4.5/xm are dominated by 
emission from the puffed-up wall at the dust sublima- 
tion radius that marks the very inner edge of the dusty 
disk. The flux from this ring is proportional to its emit- 
ting area times the blackbody function evaluated at the 
dust sublimation temperature. Determining the source 
of the variability reduces to a determination of why ei- 
ther the emitting area or temperature of the inner wall 
would change on day to week long timescales. 

Since the dust in the inner wall is already at the sub- 
limation temperature, any increase in temperature will 
cause the dust grains to evaporate. As an example, the 
timescale for a molten piece of magnesium- iron silicate 
to evaporate is given by: 



PcRc ( 27rfcT 



1/2 



(8) 



where P v is the vapor pressure ldyne cm~ 2 ), R c is 
the radius of the molten droplet, p c is the density of 
the dust grain, m is its mass and T is its temperature. 
For Twl700 K, R c «l/im, p c ~2 g cm" 3 the ey apora- 
tion time is less than a few hours (jShu et al.l lT996). The 
vapor pressure is in fact highly temperature dependent, 
inc reasing by two orders of magnitude from 1500 to 2000 
K (|Nuth fc Fergusor]|2006l ). This timescale is effectively 
instantaneous relative to our observing cadence, allowing 
us to assume that if the heating at a location on the disk 
changes it will rapidly equilibrate thermally. Assuming 
that the condensation timescale is similar to the evapo- 
ration timescale then as the heating of the disk decreases 
the inner disk will rapidly fill in smaller radii with dust 
until it reaches 1500 K. 
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Table 1 

Periodic Stars 



Star ID [3.6] period [3.6] significance [4.5] period [4.5] significance Optical Period 



Periodic stars 



LRLL 15 


15.5 


4.4 


14.8 


4.9 


not periodic 


LRLL 29 


2.2 


2.8 


2.2 


5.4 


10.8 


LRLL 36 


14.7 


4.9 


14.7 


5.54 


5.1 


LRLL 64 


8.3 


6.1 






8.1 


LRLL 69 


8.9 


6.4 


15.5 


3^6 


9.1 


LRLL 87 


12.9 


5.9 






14 


LRLL 90 


19.8 


4.2 


20.6 


4.3 


2.2 


LRLL 94 


5.6 


1.4 


5.5 


1.9 




LRLL 116 


7.6 


6.4 


4.2 


2.6 


7.0 


LRLL 125 


2.5 


2.3 


8.4 


5.6 


8.0 


LRLL 186 


8.3 


4.5 


8.3 


5.9 


not periodic 


LRLL 435 


13.9 


5.4 


13.6 


1.3 




LRLL 54361 


21. 5 a 


6.7 


21.4 a 


6.7 




Likely Periodic stars 


LRLL 41 


12.5 


3.4 


12.3 


3.6 


2.8 


LRLL 47 


4.9 


2.7 


4.9 


2.9 


4.857 


LRLL 53 


3.0 


3.2 


2.9 


3.6 


3.0 


LRLL 60A 


3.7 


4.1 


3.8 


2.2 


6.340 


LRLL 79 


1.9 


3.4 


1.9 


3.3 


1.9 


LRLL 85 


2.7 


1.4 






2.7 


LRLL 88 


5.2 


3.3 


2.5 


1.4 


5.5 


LRLL 93 


4.3 


3.5 






1.5 


LRLL 149 


2.5 


2.1 


2.5 


4.8 


2.5 


LRLL 178 


6.8 


3.4 






6.9 


LRLL 182 


2.7 


3.3 


2.6 


2.5 


2.7 


LRLL 228 


4.8 


4.0 


4.8 


3.6 


not periodic 


LRLL 276 


4.6 


3.9 


1.6 


1.8 




LRLL 40182 


2.6 


3.2 


2.6 


3.7 





Note. — The significance corresponds to the log of the probability that the peak is due to 
random fluctuations. A value of 4 corresponds to a 1 in 10 4 chance that the peak in the SigSpec 
analysis is due to noise. 

a The period presented here is derived using only the fall 2009 Spitzer data. Analysis of all 
available infrared data indicates a period closer to 25 days (Muzerolle et al. submitted) 



loose constraints on the typical size of the structural 
fluctuation. As noted earlier, the size of the [3.6] and 
[4.5] fluctuations depends on the effective temperature 
of the star, since the [3.6] and [4.5] fluxes arise partly 
from the photosphere; we account for this effect by using 
the known effective temperatures of the variable class II 
sources to simulate the observed fluctuations. For each 
star in our sample we start with blackbody emission at 
the stellar effective temperate plus a blackbody repre- 
senting the dust at 1500 K with a flux ratio of the disk 
and star emission at K band of Fdi S fe/F stor =0.4. We as- 
sume that fd and /n are normally distributed about a 
mean value with some small dispersion and draw ran- 
dom values for each star from these distributions. New 
disk fluxes at 3.6 and 4.5^tm are calculated and added to 
the stellar flux at these wavelengths. The ratios between 
the old and new 3.6 and 4.5/xm stellar+disk fluxes are 
converted to magnitudes and compared to the observed 
distribution of max-min magnitudes in each band. This 
simulation allows us to estimate appropriate values of /q 
for a given value of fd- 

As expected, fd and /n are highly anti-correlated. If 
the entire inner disk is subject to a structural perturba- 
tion (fd = 1) then fa ~ 1.25, a 25% increase in emitting 
area. This increase in emitting area could correspond to 
an increase in the radius of the inner disk, an increase in 
the scale height or some combination of the two. When 



That is, on our observing cadence the dust should al- 
ways be near 1500 K, consistent with recent observations 
of the t emperature of the inner rim as its strength fluc- 
tuates (jFlahertv et all 1201 11 I2012D . This only leaves a 
change in emitting area as the possible source of the vari- 
ability. We consider a simple scenario in which some frac- 
tion of the inner wall, fd, undergoes a fractional change 
in its emitting area, quantified as /q. The old disk emis- 
sion depends on the fiducial emitting area, while the new 
disk emission depends on both fd and /q. 

F old cx RoH B x {T d ) (9) 
F new cx ((1 - f d )R Ho + fdfnRoHo)B x (T d ) (10) 
cx ((1 - f d ) + fdfa)F a id (11) 

The parameters Rq and Hq are the fiducial radius and 
scale height of the inner disk while B\(Td) is the black- 
body function evaluated at the dust sublimation tem- 
perature. The change in disk emission depends on both 
the fractional change level of emission and the amount 
of the disk that is subject to this perturbation. By fo- 
cusing on changes in disk emission, rather than the abso- 
lute level of disk flux, we can ignore complicating factors 
such as the dust grain size and composition (which set 
the fi ducial values of Ro and H (|Dullemond fc Monnierl 
2010)). While fd and fa are highly degenerate, we can 
vary /n assuming a value for fd until the difference be- 
tween F new and F a id matches the observations, to put 
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only 1% of the inner wall undergoes a structural vari- 
ation (fd — 0.01), then /q ^25. If the radius of the 
inner wall stayed fixed, then this would imply that the 
scale height of 1% of the inner disk increases by over 
an order of magnitude. If this change in emitting area 
seems extreme, then instead a substantial portion of the 
disk must be undergoing a structural perturbation. If 
only a few tens of percent of the inner wall fluctuates 
(fd = 0.1 — 0.3) then the emitting area must change by 
a factor of 2-3. 

If 'dippers' and UXors exhibit structural perturbations 
similar to those observed here then we can use their ob- 
servations to help constrain f^ and fh. The short-lived 
extinction events in these sources are likely due to clouds 
above the midplane near the inner edge of the disk. In 
our nomenclature one could say that f^ of the disk has 
had its scale height increased by fn ■ The ratio of the du- 
ration of the extinction eve nt to the period of t he event 
corresponds to f^. AA Tau (jBouvier et al.|[2007l ) exhibits 
an extinction event that occurs every 8 days and lasts 
~1.5 days co rresponding to fd ~ 0.2 Based on the ob- 
servations of IMorales-C aldcron et al.l (|2011[ ) typical ex- 
tinction events last around 1 day with periods between 
2 and 10 days, corresponding to f<j ~ 0.1 — 0.5. 

One key distinction between the dippers and our ob- 
servations is that the dippers can be explained by a long- 
lived structural perturbation that rotates with the rest of 
the disk while our observations require a rapid change in 
the structure. A non-axisymmetric perturbation, such as 
a warp, rotating with the disk pr oduces infrared fluctu- 
ation s on the order of only 1-2% (Fl aherty fc Muzerolld 
2010), much smaller than our observations. To explain 
the variations we see, fn must represent how much ei- 
ther the radius or height of the inner disk changes on 
weekly to monthly timescales. The fact that the extinc- 
tion events in AA T au are seen in obse rvations sepa- 
rated by many years (|Bouvier et al J 120071) indicates that 
the warp causing the extinction is fairly stable. How- 
ever, n ot all of the dippers in Orion exhibit periodic be- 
havior (jMorales-Calderon et al.f2011l ) , suggesting that at 
least some of these structural perturbations are short- 
lived and may be related to the infrared variability. The 
dipper timescale then gives some idea of the size of the 
typical disk structure that might vary, although more ob- 
servations are needed to directly connect these two phe- 
nomena. 

One way to increase the emitting area is to increase the 
amount of dust in the inner disk. Increasing the surface 
density will increase the optical depth causing the r = 1 
surface, which defines the size of the emitting region, to 
occur at a larger height above the midplane. The ac- 
cretion flow can increase the surface density by dragging 
dust in and out of the disk, and it does so on a viscous 
timescale. At 0.1 AU the viscous timescale is on the order 
of thousands to tens of thousands of years, much longer 
than the weekly fluctuations that we observe, making it 
unlikely that a change in the accretion flow is causing an 
observable change in surface density. A disk wind could 
also pull material out of the disk, decreasing the surface 
density, causing the disk to appear to shrink. Near the 
base of the wind, at th e surface of the disk the s peed of 
the wind is sub-sonic (jBlandford fc Pavndfl~982T) . which 
implies a timescale to remove material that is longer than 
the thermal timescale. The thermal timescale is on the 



order of years to decades, implying that a wind could 
not remove dust fast enough to match the observed fluc- 
tuations. Given the long timescales for the accumula- 
tion/removal of dust by both accretion and a wind it 
is unlikely that we are observing a large change in dust 
mass. 

A number of other theories, besides a variable ac- 
cretion flow through the disk and a variable wind, 
can produce fluctuations in disk structure on daily t o 
weekly times c ales ( see figure 12 in lFlahertv et al.ll201lD. 
iTurner et al.l (|201Q[ ) and more recently iHirose fc Turner! 
(|2011h . find that in regions of the disk unstable to 
magneto-rotational instability (MRI) the disk magnetic 
field can become buoyant and lift out of the midplane. 
The magnetic fields drags dust along with it, caus- 
ing rapid fluctu ations in the scale height of the disk. 
iKe et all (|2012h model irradiation by X-ray flares and 
find that the the increased ionization caused by an X- 
ray flare will lead to dust being accelerated up along 
the stellar magnetic field lines away from the midplane. 
This would lead to changes in dust scale height that are 
directly correlated with large X-ray flares. If the cou- 
pling between the stellar magnetic field and the disk is 
moderate then the stellar magnetic field will periodically 
inflat e, possibly leading to the c reation of a warp in the 
disk (|Goodson fc Wingled [1999V A companion circling 
the primary star on an orbit that is inclined with re- 
spect to the disk can create a warp as it drag s material 
out of the midplane (jFragner & Nelson! l2010T. Another 
possibility, which we consider in detail below, is variable 
illumination by a hot spot and the accompanying change 
in structure as it sublimates part of the inner wall. 

5.2. Variable Illumination from a Stellar Hot Spot 

Variable illumination from a hot spot on the surface 
of the star is one possible mech anism for explaining the 
obser ved infrared variability (fMoralcs-Caldcro n et al.l 
2009). As gas flows along the stellar magnetic fields lines 
it creates a shock as it strikes the surface of the star. 
The emission from this shocked-material is strongest at 
shorter wavelengths due to its high temperature and can 
be approximated with blackb ody emission at 6000-10000 
K (|Calvet fc Gullbriridri998ft . The increased luminosity 
from this hot spot is responsible for the large irregular 
optical fluctuations seen among actively accreting stars 
(jHerbst et al.l !T994). The lack of periodicity in the opti- 
cal suggests that the distribution of spots on the stellar 
surface is asymmetric and that the spots do not live for 
multiple rotation periods. A rapid change in the covering 
fraction of these hot spots would lead to rapid variations 
in the illumination of the disk. Here we examine the 
effect of variable illumination on the disk structure to 
determine if it can plausibly explain our observations. 

We restrict ourselves to an examination of the emission 
from the directly illuminated puffed-up wall at the sub- 
limation radius that defines the inner edge of the disk, 
since the 3.6 and 4.5/mi excess flux is dominated by emis- 
sion from this region. As discussed earlier, if the inner 
wall is heated above 1500 K the dust will rapidly evap- 
orate and the wall will push outward until it reaches a 
radius defined by the new sublimation radius. This sud- 
den loss of dust over the newly illuminated region of the 
disk may have an effect on the dynamics of the system. 
While the gas dominates the mass, and hence the dy- 
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namics, the dust dominates the opacity and is responsi- 
ble for setting the disk temperature and the gas may be 
affected by this loss of opacity. Without the dust, the 
gas in this azimuthal region of the disk will cool off and 
collapse to the midplane on a thermal timescale. The 
cooling timescale is roughly the thermal energy (=nkT) 
divided by the cooling rate (=A). For the temperatures 
in the inner disk, the cooling rate is A/n e n p sa 10 -26 erg 
cm 3 s _1 . Assuming the surface density is 1000 g cm~ 2 , 
the height of the disk is 1 AU, and the disk is made en- 
tirely of hydrogen molecules, we obtain n~ 10 13 cm~ 3 . 
Assuming the ionization fraction is ^10~ 12 , to convert 
from n/i to n e gives a cooling timescale of tens of thou- 
sands of years. This is much longer than the lifetime of 
the hot spot, and as a result, the gas will not instantly 
collapse to the midplane when the dust disappears. If 
the condensation timescale is not substantially different 
from the evaporation timescale, then once the hot spot 
rotates out of view the dust will quickly fill in the disk 
back to the original sublimation radius. 

The fact that the dust can sublimate fairly quickly and 
that the gas, which drives the dynamics of the disk, re- 
sponds more slowly means that a hot spot rotating across 
the surface of the star would sweep across the inner wall 
like a lighthouse, creating a small region that has a higher 
sublimation radius than the rest of the wall. This kink 
in the inner wall (Fig |23| could lead to a change in the 
strength of the infrared emission. Given that hot spots 
are not very long lived, the size of the kink will constantly 
fluctuate on timescales of days to weeks, similar to the 
infrared variability timescales. 

The next question is whether the changes induced by 
variable illumination are large enough to produce the 
observ ed infrared fluctuations. iMorales- Caldcro n et al.l 
(2009) decided that a hot spot alone could not produce 
large enough infrared variations. They had to include a 
warp in the disk that rotated in step with the hot spot to 
create fluctuations of ~ 0.1 mag. Here we consider the 
effect of the changing sublimation radius in response to 
the hot spot to see if this structural perturbation is large 
enough to match the observations. As discussed earlier, 
the emission from the rim is proportional to its emitting 
area. 

H 

F rim CX RrirnH r i m B v {Td) = — R rim B v (Td) (12) 

If we assume that H/R stays roughly constant as the 
disk is pushed outward, then the flux from the rim 
is proportional to R 2 . The radius of the inner disk 
roughly scales as the square ro ot of the stellar luminosity 
(Dullemond & Monnicr 2010), and in turn the inner rim 
flux scales directly with the stellar luminosity. A 10,000 
K spot covering 1% of the stellar surface of a 4,000K 
star would produce stellar luminosity fluctuations close 
to 30% and a corresponding 30% increase in disk flux. 
This is roughly consistent with the size of the infrared 
fluctuations that we have observed. 

We now estimate the distribution of infrared variability 
among the cluster members in the hot spot model. The 
observed changes in [3.6] and [4.5] will depend on the size 
of the hot spot and on the temperature of the star. A hot 
spot covering 1% of a M dwarf will have a more significant 
effect on the luminosity striking the disk than if the same 



hot spot covered 1% of a G star because of the larger tem- 
perature contrast. Also, an M star has a stellar effective 
temperature closer to the dust temperature, leading to 
smaller changes in [3.6] and [4.5] for a given change in 
disk flux. Our estimate is based on the known distribu- 
tion of effective temperatures for the cluster members, 
while also making some simple assumptions about the 
typical size of a hot spot, its temperature, how much of 
the disk is heated by the hot spot and the strength of 
the excess. We start with the effective temperatures of 
the variable stars with — 1.8 < ajRAC < 0. For the disk 
we use a 1500 K blackbody that has Fdisk/F star =0-4 at 
K-band, and sum the stellar and disk flux at 3.6 and 
4.5/im to determine the fiducial system flux. We ran- 
domly draw a hot spot size from a normal distribution 
with a given mean and standard deviation and add a spot 
of this size with T=8000 K to the stellar flux. We assume 
that only one-third of the disk is illuminated by the hot 
spot (fd=l/3); the flux from the other two-thirds of the 
disk stays constant. The illuminated region of the disk 
has its flux increased by the same factor as the change in 
stellar luminosity for the given spot size and stellar effec- 
tive temperature (fn = AL*/L*). We then add this new 
dust flux to the new stellar flux at 3.6 and 4.5/mi and 
calculate the ratio of this new flux to the fiducial flux. 
We convert the flux ratios at [3.6] and [4.5] into mag- 
nitude differences and compare the distribution of these 
differences to the observations. We also compare our 
model with observed I-band fluctuation s in the irregu- 
lar variables in IC348 (jCohen et al.ll20 04). Including the 
optical data, while not available for all of the infrared 
variables, helps constrain the size of the star spot inde- 
pendent of our assumptions about the disk. Matching 
to the observed fluctuations at all wavelengths allows us 
to determine if a consistent set of hot spot sizes can re- 
produce the observed distribution of fluctuations. Given 
that the size of each hot spot is random, many simula- 
tions are repeated in order to determine how often the 
distribution of magnitude fluctuations matches the ob- 
servations. The uncertainty in the assumptions (ie. hot 
spot temperature, how much of the disk is heated by the 
hot spot, the exact response of disk flux to a change in 
stellar flux) prevent us from using this simulation to ac- 
curately constrain the hot spot sizes, but it can be used 
to check for plausible solutions. 

We find that hot spots that are 8000 K, with a mean 
size of 3±1% of the stellar surface are able to reproduce 
the observed fluctuations in [3.6] and [4.5] as well as I 
band. We can also fit the observations with a 10,000K 
hot spot covering on average 2.3±0.7% of the surface and 
heating 15% of the disk (Fig I22|) . Both of these scenar- 
ios produce infrared fluctuations of ^0.2 mag as well as 
optical fluctuations of 0.5-1.5 mag, similar to the obser- 
vations. This suggests that heating from a hot spot is 
a reasonable model for t he fluctuations. The structura l 
perturbations required bv lMorales-Calderon et all (2009) 
can naturally be explained by the increase in disk radius 
as dust rapidly sublimates in response to the increased 
illumination. We do assume that H/R stays constant as 
the inner edge of the disk is pushed outward, but this as- 
sumption can be relaxed if f^ is allowed to vary. If H/R 
drops by a factor of 2, then f^ needs to increase by a factor 
of 2 from 1/3 (15%) to 2/3 (30%) for an 8,000 (10,000) K 
hot spot to fit the observations. The change in H/R and 
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frf are highly degenerate, preventing us from accurately 
constraining either of them. Once the material behind 
the inner wall becomes directly illuminated by the star 
it will heat up, but its gas pressure scale height will not 
change substantially because the thermal timescale is on 
the order of years and decades. The observed height will 
simply be the fiducial height. The region directly behind 
the inner wall lies in the shadow of the wall and will 
be heated by radial diffusion which is difficult to model 
(jDullemond fc Monnierl 120101 ). Our assumption of con- 
stant H/R is most likely not accurate, but is simply used 
as a first order approximation of reality. The fact that 
these first order results are reasonable suggest that this 
model warrants further investigation, including a better 
treatment of the height of the disk. 

5.2.1. How well does the hot spot model fit the trends we see 

in our data? 

We have found that infrared variability is strongly cor- 
related with large irregular optical fluctuations, position 
within the cluster and X-ray luminosity. Here we ex- 
amine how those correlations could be explained with 
illumination of the disk by variable hot spots. 

Optical Variability: Large irregular optical fluctuations 
are directly attributed to the presence of hot spots on the 
stellar surface. The strong influence of hot spots on disk 
structure is consistent with our finding that almost 90% 
of the class II sources with large irregular optical fluctu- 
ations are also variable in the infrared while only 30% of 
the class II sources without irregular optical fluctuations 
are variable in the infrared. While the optical data were 
taken six years before our infrared data, it is plausible 
that the size of hot spot fluctuations does not change over 
the course of a few years. While the short-lifetime of hot 
spots may make it difficult to search for a direct correla- 
tion between optical and infrared fluctuations, contem- 
poraneous observations may show a correlation between 
the size of the infrared and optical fluctuations. 

The South- West Region: Along the southwestern edge 
of the main cluster is a region of higher density gas and 
dust, in which many of the mos t deeply embe dded mem- 
bers of this cluster are located (lHerbstll2008ft. The mam 
concentratio n of sub- mm sources (| Jorgenson et al.f 2008) , 
H2 outflows (Eisl offel et al.l 12003ft and class I and class 
objects (|Muench et al.l 12007ft are in this southwestern 
ridge , which has a total mass of - 10M Q (jTafalla et al.l 
2006). The presence of stars still embedded in massive 
infalling envelopes, and the youn g dynamical age of the 
outflows (|Walawander et aLll2006l ). indicate that this is a 
region of very recent star formation. The relative youth 
of this site may be the most relevant feature for the vari- 
ability. IC 348 is not home to any massive O stars that 
could be producing strong photo-ionizing radiation and 
we do not expect that the ambient radiation field will 
vary strongly across the cluster. While infrared classes 
are broadly associated with the age of a source, it is pos- 
sible that within a single SED class there is a spread in 
ages. The class II sources in the south-western ridge may 
be younger that the class II sources in the more exposed 
cluster center. Their relative youth may be associated 
with higher av erage accretion rates. A ccretion rates de- 
crease as t~ 2 (Hartma nn et al.l 11998ft and from 1 to 3 
Myr, a plausible age spread between regions in IC 348, 
the average accretion rate will decrease by a factor of six. 



iCalvet fc Gullbrinel ([1998ft find that the excess from the 
hot spots around stars with different accretion rate can 
be explained by simply changing the covering area of the 
hot spot, with little change in the density and temper- 
ature within the spot. Stars with higher accretion rates 
have more of their surface covered by hot spots. Simu- 
lations find that the flow of gas onto the star is highly 
variable w ith more severe instabilities for higher accre- 
tion rates (jKulkarni fc Romanovall2008l ) and it is possible 
that a larger average hot spot size, due to a higher accre- 
tion rates, would be associated with larger fluctuations 
in their coverage. In this way, the higher accretion rate 
of the younger sources in the southwestern ridge may be 
associated with larger more variable spots that in turn 
lead to larger disk structural fluctuations. 

The importance of age would imply that class I sources 
(aiRAC > 0) should be more variable that class II sources 
(—1.8 < aiBAC < 0) since class I sources represent an 
earlier stage of stellar evolution. In reality we do not see a 
trend of increasing f var from class II to class I sources but 
this may partly be due to the fact that amAC is an im- 
perfect measure of evolutionary stage. Some nearly edge- 
on class II sources can appear to have a sharply rising 
infrared SED, while nearly pole-on class I sources have a 
spect ral shape similar to class II sources ([Robitaill e et al.l 
[2007ft . This mixing of evolutionary stages will dilute any 
trend between f var and SED shape among the youngest 
sources. It is also possible that our limitation of only con- 
sidering sources with [3.6] < 14 in our statistical analysis 
preferentially selects older class I sources with smaller 
envelopes, and hence less extinction of the central disk 
emitting at 3/im, further diluting our results. 

X-ray Luminosity: The connection between X-ray 
emission and accretion has been difficult to establish, 
producing conflicting results. Actively accreting stars are 
known to be on average less X-ray lumi nous than stars 
that have no strong ongoing accretion (Fcigclso n et al.l 
12007ft . In some stars the soft X-ray component has 
been associated with low temperature, high-density ma- 
terial in the accretion shock, while the hard X-ray 
flux arises from hotter, low-density g as in the corona 
([Argiroffi et al.l 12011ft . iStassun et aTT (|2007ft find that 
optical variability is strongly correlated with X-ray lumi- 
nosity, which they suggest may be due to larger magnetic 
topology giving rise to both stronger X-ray emission and 
larger spots. In terms of our analysis, this implies that 
stars that are more X-ray luminous have larger spots, 
giving rise to larger fluctuations in disk illumination and 
disk structure. The lack of correlation among the class 
III sources is simply due to the fact that these stars do 
not have any dust that can reflect the varying spot struc- 
ture. 

5.2.2. How well does the hot spot model fit other literature 

results ? 

In our analysis, we restricted ourselves to the puffed- 
up inner rim, but an extended radial region of the disk 
will be heated to a higher than normal temperature by 
the hot spot. Longer wavelength flux arising from large 
radii would rise in concert with the short-wavelength 
flux obse rved here, as has b een seen in infrared spec- 
troscopy (Kospal ct al. 2012) and photometry (Megeath 
et al. accepted) of class II sources. Both of these studies 
found that for the majority of class II sources the increase 
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Figure 22. Comparison of magnitude fluctuations derived from out 10,000 K hot spot model (black solid line) and the observed fluctuations 
for the variable class II sources (red dotted line). We show, from top to bottom, the max-min [3.6], [4.5] and I band magnitude. F o r the I 
band data the number of counts have been normalized since we use only the 15 stars observed as irregular variables by Cohen et al. (2004), 
but include the entire simulated sample (~50 stars). Dashed lines include the sources that would lie below our infrared detection limit. 
Adding 10,000 K hot spots with an average size of 2.5% of the stellar surface, and accounting for the change in inner radius due to the 
additional heating from these hot spots, is able to reproduce the observed infrared fluctuations. Similar results are found when using an 
8,000 K hot spot with an average size of 3% of the stellar surface. 



in 3-10^m flux was relatively constant with wavelength. 
These wavelengths are sensitive to the surface layers of 
the disk, which respond almost instantaneously t o an 
increase in illumination ([Chiang fc GoldreicldlT997r i. al- 
lowing for the disk flux to increase as the illumination 
increases. 

This wavelength dependence contrasts with many tran- 
sition disks that have infrared fluctuations in which 
the short- wavelength (A < 8/xm) flux decreases while 
the long- w avelength (A > 8/zm ) flux increases , and 
vice-versa dMuzerolle et all l2009t lEspaillat et~"aH l2011t 
iFlahertv et al.l 120121) . This 'seesaw' behavior has been 
explained by changing the scale height of the puffcd- 
up inner wall; as the scale height increases the short- 
wavelength flux emitted by this part of the disk increases 
while the long-wavelength flux decreases as the outer 
disk becomes more shadowed by the inner disk. This 
model has been successful at matching the observations 
but it may not be the only explanation. The observa- 
tions themselves do not necessarily require a change in 
the scale height; they only require that the covering frac- 
tion of the inner disk as seen from the outer disk changes 
with time. In an axisymmetric model the only way to 
accomplish this is to change the scale height of the in- 
ner disk. Another possibility is that the scale height 
stays constant, but the azimuthal extent of the inner wall 
changes. Instead of being a solid ring wrapped around 
the star, the wall may be a crescent whose size changes 
with time. In the context of our hot spot model, this 
could happen as an accretion spot heats and sublimates 
a portion of the inner wall. If the wall is radially thin and 
there is no dust behi nd it, a geometry tha t is consistent 
with the SED shapes (lEspa iflat et aljfeoiOl ). then instead 



of putting a dent in the wall the increased luminosity 
from the hot spot will punch a hole through it (Fig |23| . 
The outer disk will become directly illuminated, increas- 
ing the long wavelength flux while the short wavelength 
flux decreases because the emitting area of the inner wall 
has decreased. In this way illumination from variable ac- 
cretion spots can explain the wavelength dependence of 
both class II and transition disk sources. 

Observations of the accretion luminosity can put lim- 
its on the application of this model. We have previously 
studied a handful of evolved disks using near-infrared 
hydrogen recombinatio n lines whose flux is propor tional 
to the accretion rate (|Flahertv et al.l 1201 ll 1201 2D . We 
found that the accretion rate was highly variable, some- 
times dropping below our detection limit, but in general 
the accretion luminosity was small compared to the stel- 
lar luminosity (L acc /L* <.06). This would imply that 
the change in luminosity due to the appearance of a hot 
spot would be small, creating a small change in inner 
wall radius. For a 6% increase in luminosity the radius 
only increases by 6%. For our model of an accretion 
spot blowing a hole in the inner wall, the radial depth 
of the wall must be very small (AR/R < .06). SED 
modeling of pre-transition disks often assumes that the 
inner wall is infinitesimally th in and is still able to fit 
the SED (|Espaillat et al.l 12010). While the presence of 
dust directly behind the puffed-up inner wall is difficult 
to constrain because it is shadowed by the wall and ap- 
pears very faint, a very narrow wall is consistent with 
the observations. In these previous studies, a direct cor- 
relation between infrared and accretion rate fluctuations 
was not seen, but this may be due to the complex ge- 
ometry of the accretion flow, in combination with our 
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Figure 23. Schematic diagrams illustrating how the presence of a hot spot can lead to large infrared fluctuations. When a hot spot 
appears on the surface of the star it heats up the disk, causing the dust to rapidly sublimate, effectively pushing out the inner edge of the 
disk (left and middle left). This increase in emitting area causes an increase in infrared flux. In the case of a transition disk, in which 
there is very little dust behind the wall, instead of creating a dent the hot spot punches a hole in the inner disk (middle right and right). 
This decreases the short-wavelength flux, while also directly illuminating the outer disk, which increases the long-wavelength flux. Given 
the small observed accretion luminosities among evolved disks, the inner wall must be very thin, but otherwise it is consistent with the 
'seesaw' behavior that is often observed in transition disks. 

particular line of sight. A hot spot on the far side of 
the star would produce a large change in infrared excess 
while not producing a significant change in the line flux. 
More detailed models, and continued observations, may 
help to resolve this potential conflict. 

One possible exce ption to our model is LRLL 2. 
lEspaillat et all ([2012) were able to model its SED with- 
out a large gap. This star still exhibits a 'seesaw' behav- 
ior in its SED, which is difficult to explain with hot spot 
illumination without a narrow wall and a gap. It is pos- 
sible that there is in fact a gap in the di sk, but it is too 
small to strongly effect the SED shape. lAndrews et al.l 
(2011) find that some sources with a gap based on re- 
solved sub-mm observations do not show the traditional 
SED shape of a pre-transition disk. It is possible that 
LRLL 2 represents another example of a system that has 
a radial gap that does not announce itself in the SED. 



5.3. Class III Variables 

For Class III sources, which by definition have little to 
no dust emission, alternative methods of producing in- 
frared variability must be considered. Out of our entire 
sample of 108 class III sources, we find 21 possible vari- 
able stars. Four sources (LRLL 154, 286, 298, 329) have 
only a single discrepant data point and are unlikely to be 
real variables. Four sources (LRLL 64, 69, 87, 125) are 
periodic in the optical and infrared, with similar periods 
in both. Periodic fluctuations in the optical are most 
likely driven by cool spots rotating across the surface of 
the star, and we are likely seeing the effects of these spots 
in the infrared. One source, LRLL 314, is periodic in the 
optical but was not detected as periodic in the infrared, 
possibly because of the weakness of the signal. Cool spots 
should produce smaller fluctuations in the infrared than 
in the optical, and it may be difficult to pick out period- 
icity given the small size of the fluctuations. One source, 
LRLL 29, has an infrared period of 2.2 days, but has an 
optical period of 1 0.8 days. While so me fluctuations in 
the period are seen (|Cohen et al.ll2004T ). these changes are 
too small to account for the difference reported here. It is 
possible that period aliasing has led to different measures 
of the period in the optical and infrared. The remain- 
ing ten sources (LRLL 1, 4, 47, 48, 53, 179, 175, 259a, 
277, 1939; see Table for positions and other ancillary 
information) have no evidence for an optical or infrared 
period. The fluctuations tend to be weak, with min to 
max differences less than 0.2 magnitudes. None of them 



are detected at [24], ruling out a cold massive disk that 
does not produce any excess over the 3 to 8/xm range. 
LRLL 1,4 53 show brief epochs with a small handful of 
data points fainter than the rest of the sample, simila r 
to the 'dippers' seen by Morales- Caldcr on et al.l (|2011f) . 
These could be extinction events caused by small clumps 
of dust, which may not be massive enough to create a 
measurable excess. Another possibility is that there ex- 
ists an unseen companion with a circumstellar disk that 
is variable. An M6 star with a disk whose K-band veil- 
ing was close to zero but occasionally increases by a few 
tenths, would produce [3.6] and [4.5] fluctuations of 0.05 
to 0.2, similar to what is observed. While the presence of 
this dust would create a small excess, increasing amAC 
by 0.15 at the most above the single star value, the in- 
frared SED may still be steep enough to be counted as 
a diskless source. More detailed followup is needed to 
confirm the variability as well as determine its origin. 

6. CONCLUSION 

We have used 38 epochs of [3.6] and [4.5] photometry 
spread over 40 days and of over 200 low-mass members 
of the IC 348 cluster to study the structure of proto- 
planetary disks. We have found that many of the stars 
are variable in the infrared and that this variability is 
strongest and most frequent among the stars with the 
strongest infrared excesses. Fluctuations of a few tenths 
of a magnitude on weekly timescales are often seen from 
these young systems. It also appears that the systems 
in the south-west ridge, which are most likely younger 
than the rest of the cluster, are more likely to be vari- 
able. Those sources with larger X-ray luminosity also 
seem to be more variable, as well as stars with large ir- 
regular optical fluctuations. Relatively few of the stars 
show periodic fluctuations, most of which are related to 
stellar fluctuations from cold spots rather than changes 
in disk structure. 

The frequency and size of the variability indicate that 
the disk is not a static entity. One likely possibility is 
that there are large structural changes in the disk, due 
to dust and gas being moved on a dynamical timescale, 
which at the inner edge of the disk is on the order of a 
week for many of these stars. We also find that dust can 
quickly sublimate when heated above 1500 K, which may 
result in large changes in the structure of the disk. Both 
the movement of material, and the rapid sublimation of 
dust, will be manifest as large structural changes in the 
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emitting area of the inner disk. This increase in emit- 
ting area must be substantial in order to produce the 
observed fluctuations; increasing the scale height of 30% 
of the disk by a factor of 2 could reproduce the observa- 
tions. A broad range of scale height increases can also 
explain the measurements if appropriate adjustments are 
made in the amount of the disk undergoing a perturba- 
tion. Multiple theories for structural perturbations exist 
in the literature, including buoyant disk magnetic fields 
lifting dust out of the midplane, X-ray flares accelerating 
dust out of the midplane, and perturbations by the stel- 
lar magnetic field. Our observations suggest that these 
theories must be capable of structural changes on large 
physical scale. We explore the effect of variable illumi- 
nation by a hot spot in additional detail. A 10,000 K 
hot spot covering on average 2.5% of the stellar surface 
and illuminating 15% of the inner wall can reproduce 
the size of the fluctuations and is consistent with the ob- 
served timescale. Systems that are younger and more X- 
ray luminous are more likely to have variable accretion, 
which would explain the ancillary trends that have been 
seen. While this model is dependent on many simplifying 
assumptions, it does indicate that variable illumination 
warrants further consideration. 

This study of infrared variability, as well as others, 
points to a new picture of the disk. It is no longer a 
static system that slowly evolves, but it is a bubbling, 
boiling, wrinkled, dented, warped mass of gas and dust. 
These asymmetric structures cover a substantial portion 
of the disk and appear/disappear over the course of a 
month. While their source is still uncertain, it is clear 
that they are common and large. The structures are too 
small to be resolved with current instruments, but with 
infrared variability we have learned about their presence. 
Future resolved observations should frequently find non- 
axisymmetric structure that may have a large impact on 
the growth and migration of planetesimals within the 
disk. 
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Figure 24. [3.6] (left) and [4.5] (right) light curves for cluster members. Error bars on individual points have been included, as well 
as the error bars including the systematic uncertainty (vertical line, description of the derivation in the text). Only data points with an 
uncertainty of less than 0.5 mag are shown. Magnitudes have been normalized to the mean magnitudes, which are marked in the y-axis 
label. We show one set of light curves here, while the rest are available as an online figure. 
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Table 2 

Cluster Member Properties 



LRLL 


RA 


Dec 


T e// 


L* 


Ha EW a 


Optical 




Lx c 


[3.6] 




atiRAC c 


# 


(J2000) 


(J2000) 




(A) 


Period b 




(10 28 


excess 








(degrees) 


(degrees) 


(K) 


(L©) 


(days) 




erg/s) 








1 


56.1425 


32.1629 


15400 


605.000 


-10.000 


-10.000 


1.910 


-10.000 


0.010 


-2.76 


-2.80 


2 


56.1474 


32.1679 


8970 


137.000 


-10.000 


0.000 


3.830 


1053.000 


0.003 


-1.41 


-1.47 


3 


56.2110 


32.3185 


9520 


135.000 


-10.000 


6.097 


4.470 


-10.000 


-0.012 


-2.77 


-2.85 


1 


56.1300 


32.1061 


7200 


34.000 


-10.000 


0.000 


1.770 


14.900 


-0.004 


-2.69 


-2.72 


5 


56.1084 


32.0751 


5520 


9.900 


5.000 


6.536 


5.320 


385.000 


0.188 


-1.14 


-1.24 


6 


56.1539 


32.1126 


5830 


17.000 


-10.000 


1.659 


3.900 


1961.000 


-0.007 


-2.02 


-2.08 


7 


56.0353 


32.1213 


9520 


35.000 


-10.000 


-10.000 


0.920 


-10.000 


-0.002 


-2.70 


-2.72 


8 


56.0381 


32.1193 


8970 


28.000 


-10.000 


-10.000 


0.960 


-10.000 


-0.001 


-2.38 


-2.40 


9 


56.1632 


32.1551 


5520 


11.000 


-10.000 


2.600 


5.530 


573.000 


-0.015 


-2.76 


-2.86 


10 


56.1028 


32.1708 


6890 


13.000 


-10.000 


0.000 


1.910 


13.200 


-0.004 


-2.86 


-2.89 


Note. — The 


full version of this table is 


available 


in machine-readable format. Here 


we show a small portion for 


illustration. The 


full tabic includes columns for LRLL RA, Dec, T e //, L*, HqEW, Optical Period, Ay 


, hx , [3.6] excess, a. 


&IRAC, 


HH 211 distance, 



Var? and Referenc e. Values of -10 or -100 correspond to a non-detection. We do not include sources 8024, 8042, 8078, 9042, 9060, 30190, 
30191, 30192 from lLada et all 11200611 because they are duplicates of sources 24A, 42A, 78A, 42B, 60A, 92,1,9 respectively 
a Here we actually report the negative of the Ha EW, with positive numbers corresponding to emission. Ha EW is not reported when it 
is in absorption, and a value of -10 corresponds to a non-detection. 

Optical period is taken from lCohen et al.l d2004D :IC ieza & Ba libcr (200l|) and references therein. Values of correspond to stars that 
were detected in these surveys, but no period structure was detected, while values of -10 refer to stars that were not reported as detections 
in these surveys. 

c X-ray luminosity taken from Prcibisch & Zinncckcr (2002) 

" The slope of the observed infrared SED, defined as AFa oc X a over the IRAC (3-8/mi) wavelengths. 

° The slope of the dereddened infrared SED. Unless otherwise stated, this is the value of the infrared SED that is used in the text. 

^ Var?: y — Variable, n — not variable, yl — Variable, but only detected at [3.6], nl — not variable, but only detected at [3.6], p = 

periodic, b — a nearby source prevents us from accurately assessing variability, ... — source was not detected 

g Ref erence key : luh03 - ILuhman et al.| 1)20031) , mue07 - IMuench et al.| | | 20071V lad06 - ILada et"al| < |2006D . muz06 — IMuzerolle et al.| 
d2006fl . Iuh05 = ILuhman et all d2005l) 



Table 3 

[3.6] photometry 



RA (J2000) 
(Degrees) 


Dec (J2000) 
(Degrees) 


55106.9 
([3.6]) 


([3.6] unc) 


55107.3 
([3.6]) 


([3.6] unc) 


55107.8 
([3.6]) 


([3.6] unc) 


55108.1 
([3.6]) 


([3.6] unc) 


56.1425 
56.1474 
56.2110 


32.1629 
32.1679 
32.3185 


6.510 
7.015 
7.575 


0.0021 
0.0018 
0.0019 


6.516 
7.009 
7.601 


0.0017 
0.0016 
0.0021 


6.512 
7.012 
7.587 


0.0021 
0.0021 
0.0019 


6.566 
6.993 
7.564 


0.0020 
0.0019 
0.0020 


Note. — The : 


full version of this 


table is 


available in 


machine-readable format. 


Here we 


show a small portion for 


illustration. 



Column headings refer to the MJD at which the data was taken. Magnitudes less than zero, with uncertainties of 10, correspond to 
a non-detection. 



Table 4 

[4.5] photometry 



RA (J2000) 


Dec (J2000) 


55106.9 




55107.3 




55107.8 




55108.1 




(Degrees) 


(Degrees) 


([4-5]) 


([4.5] unc) 


([4.5]) 


([4.5] unc) 


([4.5]) 


([4.5] unc) 


([4-5]) 


([4.5] unc) 


56.1425 


32.1629 


6.489 


0.0022 


6.507 


0.0021 


6.504 


0.0023 


6.578 


0.0022 


56.1474 


32.1679 


6.862 


0.0022 


6.863 


0.0022 


6.833 


0.0020 


6.870 


0.0022 


56.2110 


32.3185 


-100.000 


10.0000 


-100.000 


10.0000 


-100.000 


10.0000 


-100.000 


10.0000 


Note. — The 


full version of this table is 


available in 


machine-readable format. 


Here we 


show a small portion for 


illustration. 



Column headings refer to the MJD at which the data was taken. Magnitudes less than zero, with uncertainties of 10, correspond to a 
non-detection. 



